/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             *
 *              GNU Lesser General Public Licence (LGPL) version 3,             *
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
#include "HiMusicFittingTask.h"
#include "HiMusicMappedData.h"
#include "HiMusicCalibratedData.h"
#include "HiMusicWaveForm.h"
#include "HiMusicWaveEvent.h"
#include "FairTask.h"
#include "FairRootManager.h"
#include "TROOT.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TLine.h"
#include "TCanvas.h"
#include "TLegend.h"
#include <vector>

HiMusicFittingTask::HiMusicFittingTask() : FairTask("HiMusicFittingTask")
{
  LOG(debug) << "Default Constructor of HiMusicFittingTask";
}

HiMusicFittingTask::~HiMusicFittingTask()
{
  LOG(debug) << "Destructor of HiMusicFittingTask";
}

InitStatus HiMusicFittingTask::Init()
{
  LOG(debug) << "Initilization of HiMusicFittingTask";

  FairRootManager *ioman = FairRootManager::Instance();
  if (!ioman)
    return kERROR;

  fHiMusicMappedArray = (TClonesArray *)ioman->GetObject("HiMusicMappedData");
  if (!fHiMusicMappedArray)
  {
    LOG(warn) << "-W- HiMusicFittingTask::Init: No Data array!";
    return kERROR;
  }

  fMusicFittingWaveDataArray = new TClonesArray("HiMusicCalibratedData");
  ioman->Register("HiMusicCalibratedData", "Calibrating", fMusicFittingWaveDataArray, kTRUE);

  return kSUCCESS;
}

InitStatus HiMusicFittingTask::ReInit()
{
  LOG(debug) << "Initilization of HiMusicFittingTask";
  return kSUCCESS;
}

void HiMusicFittingTask::Exec(Option_t * /*option*/)
{
  LOG(debug) << "Exec of HiMusicFittingTask";

  auto *mappedData = (HiMusicMappedData *)fHiMusicMappedArray->At(0);
  ts_pre = ts;
  ts = mappedData->GetTimeSliceID();
  TClonesArray *mappedDataChannels = mappedData->GetChannels();

  int nChannels = mappedData->GetNumOfChannels();
  if (nChannels > 9)
  {
    LOG(warn) << "Number of channels exceeds 9! Fitting method only supports up to 9 channels.";
    nChannels = 9;
  }

  fMusicFittingWaveDataArray->Clear("C");
  auto *calibrated = new ((*fMusicFittingWaveDataArray)[0]) HiMusicCalibratedData();
  calibrated->SetEventID(mappedData->GetEventID());
  calibrated->SetTimeStamp(ts);
  calibrated->SetTimeStamp_pre(ts_pre);
  calibrated->SetDetID(mappedData->GetDetID());
  for (int chId = 0; chId < nChannels; ++chId)
  {
    HiMusicWaveForm *wave = (HiMusicWaveForm *)mappedDataChannels->At(chId);
    TClonesArray *fittedPeaks = FitPeaks(*wave);
    calibrated->SetChannel(chId, fittedPeaks);
  }
}

void HiMusicFittingTask::SetParContainers()
{
    // LOG(info) << "SetParContainers of NewTask";
}
void HiMusicFittingTask::Finish() {}

TClonesArray *HiMusicFittingTask::FitPeaks(HiMusicWaveForm &wave)
{

  double tStep = wave.GetTStep();
  TClonesArray *fitwaves = new TClonesArray("HiMusicWaveEvent");
  int nPoints = wave.GetNumOfPoints();
  // ----------------------------------------------------------
  // 2. Run FindPeaks to get baseline and peak positions
  // ----------------------------------------------------------
  std::vector<int> peaks;
  int nPeaks;
  float baseline;
  float noise;
  FindPeaks(wave, peaks, nPeaks, baseline, noise); // baseline and noise computed inside
  if (nPeaks == 0)
    return fitwaves;
  const float *raw = wave.GetData();

  // ----------------------------------------------------------
  // 3. Subtract baseline
  // ----------------------------------------------------------
  for (int i = 0; i < nPoints; i++)
    wave.SetData(i, raw[i] - baseline);

  // ----------------------------------------------------------
  // 4. Prepare graphics
  // ----------------------------------------------------------
  // TCanvas *can = new TCanvas("PlotWave","PlotWave",900,900);
  // can->Divide(3,3);
  // can->cd(1); wave.Plot();
  TCanvas *c = (TCanvas *)gROOT->FindObject("c");
  if (!c)
    c = new TCanvas("c", Form("%s fit", fitType.Data()), 1000, 600);
  TGraph *gOrig = new TGraph();
  gOrig->SetTitle("Original waveform");
  gOrig->SetLineColor(kBlack);
  TGraph *gResidul = new TGraph();
  gResidul->SetTitle("Residual waveform");
  gResidul->SetLineColor(kBlack);
  TGraph *gSum = new TGraph();
  gSum->SetTitle("Sum waveform");
  gSum->SetLineColor(kRed);

  for (int i = 0; i < nPoints; i++)
  {
    gOrig->SetPoint(i, i * tStep, wave.GetData(i));
    gResidul->SetPoint(i, i * tStep, wave.GetData(i));
  }

  // ----------------------------------------------------------
  // 5. Fit loop
  // ----------------------------------------------------------
  // ---------- 0. 评估第一峰之前的信号影响 ----------
  const float SLOPE_TAIL = -0.005;

  std::vector<HiMusicWaveForm> waves;
  float slope = wave.GetPedSlope(100);
  if (slope < SLOPE_TAIL && peaks[0] > 100) /// 拟合前一个信号的尾巴
  {
    LOG(debug1) << "slope<" << SLOPE_TAIL << ", 拟合前一个信号的尾巴";
    double A0 = wave.GetData(0);
    double tp = 0;
    double tao = fittingTao;
    double exp1_par[3] = {A0, tp, tao}; /// single exp_noped
    TF1 *ftail = new TF1("FitSingleExp", HiMusicFittingTask::singleExp, 0, 1000 * tStep, 3);
    ftail->SetParameters(exp1_par);
    ftail->SetParNames("pulse height", "tp", "tao");
    ftail->SetLineColor(kGreen);
    gOrig->Fit(ftail, "QNR");
    HiMusicWaveForm awave = CreateWave(ftail, nPoints, tStep);
    waves.push_back(awave);
    wave.Substract(awave); /// 扣除尾巴信号的影响
    // can->cd(2); wave.Plot();awave.Plot(2,"same");
    for (int i = 0; i < nPoints; i++)
    {
      gResidul->SetPoint(i, i * tStep, wave.GetData(i));
    }
  }

  // ---------- 1. 拟合----------
  // ----------  预计算所有峰-谷时间 ----------
  std::vector<double> valleyTime(nPeaks, 0.0);
  for (int ip = 1; ip < nPeaks; ++ip)
  {
    int prevPeakBin = peaks[ip - 1];
    int thisPeakBin = peaks[ip];
    int valleyBin = prevPeakBin;
    float valleyVal = wave.GetData(prevPeakBin);
    for (int j = prevPeakBin + 1; j <= thisPeakBin; ++j)
    {
      if (wave.GetData(j) < valleyVal)
      {
        valleyVal = wave.GetData(j);
        valleyBin = j;
      }
    }
    valleyTime[ip] = valleyBin * tStep; // 第 ip 个峰与其前一个峰之间的谷底时间
  }
  std::vector<double> amps, starts, pt;
  ;
  const double tauRise = 20.0;        // ns
  const double tauDecay = fittingTao; // ns

  std::vector<TF1 *> fits;
  for (int ip = 0; ip < nPeaks; ++ip)
  {
    int p = peaks[ip];
    double tPeak = p * tStep;

    // define fit window
    // 拟合下限: 峰前1us内10倍 噪声阈值或峰前 10τRise
    double tMin = tPeak - 10 * tauRise;
    double tMin_th = tMin;
    int tlow;
    if (p < 40)
      tlow = 0;
    else
      tlow = p - 40;
    for (int j = tlow; j < p; j++)
    {
      if (wave.GetData(j) > 10 * noise)
      {
        tMin_th = j * tStep;
        break;
      }
    }
    tMin = TMath::Max(tMin, 0.0);
    tMin = TMath::Min(tMin, tMin_th);

    //  拟合上限：下一个峰谷 or 峰后 10τDecay
    double tMax = tPeak + 10 * tauDecay;
    if (ip + 1 < nPeaks)
      tMax = TMath::Min(tMax, valleyTime[ip + 1]);

    // create TF1
    TF1 *f;
    if (fitType.EqualTo("biExp"))
    {
      f = new TF1("FitBiExp", HiMusicFittingTask::biExp, tMin_th, tMax, 4);

      f->SetParNames("Amplitude", "t0", "tauRise", "tauDecay");
      f->SetParameter(0, wave.GetData(p)); // amplitude
      f->SetParameter(1, tPeak - tauRise); // t0 (before peak)
      // f->SetParameter(1, tMin);       // t0 (before peak)
      f->SetParameter(2, tauRise);  // τrise fixed
      f->SetParameter(3, tauDecay); // τdecay
      f->SetParLimits(0, 0.0, 16384);
      f->SetParLimits(1, tPeak - 1000, tPeak);
      f->SetParLimits(3, 200.0, 2000.0);
    }
    else if (fitType.EqualTo("lineExp"))
    {
      f = new TF1("FitLineExp", HiMusicFittingTask::lineExp, tMin, tMax, 4);

      f->SetParNames("Amplitude", "t0", "tp", "tauDecay");
      f->SetParameter(0, wave.GetData(p)); // amplitude
      // f->SetParameter(1, tPeak - 0.8*tauRise);       // t0 (before peak)
      f->SetParameter(1, tMin_th);  // t0 (before peak)
      f->SetParameter(2, tPeak);    // tp
      f->SetParameter(3, tauDecay); // τdecay
                                    // f->SetParLimits(0, 0.0, 16384);
                                    // f->SetParLimits(1, tPeak-1000, tPeak);
                                    // f->SetParLimits(3, 200.0, 2000.0);
    }
    else
    {
      LOG(error) << "fit type error !!! Only support biExp and lineExp.";
      return 0;
    }
    gResidul->Fit(f, "QNR", "", tMin, tMax);
    fits.push_back(f);

    double amp = f->GetParameter(0);
    double t0 = f->GetParameter(1);
    double tp = f->GetParameter(2);
    double tao = f->GetParameter(3);

    amps.push_back(amp);
    starts.push_back(t0);
    pt.push_back(f->GetMaximumX());
    HiMusicWaveEvent wavedata;
    wavedata.SetPulseHeight(amp);
    wavedata.SetArrivalTime(t0);
    wavedata.SetRisingTime(tp - t0);
    wavedata.SetDecayTime(tao);
    wavedata.SetPedestal(baseline);
    wavedata.SetChi2OvNdf(f->GetChisquare() / f->GetNDF());

    new ((*fitwaves)[ip]) HiMusicWaveEvent(wavedata);

    // subtract fit from working waveform
    HiMusicWaveForm awave = CreateWave(f, nPoints, tStep);
    waves.push_back(awave);
    wave.Substract(awave);
    // if(ip<7){ can->cd(ip+3); wave.Plot();}
    for (int i = 0; i < nPoints; i++)
    {
      gResidul->SetPoint(i, i * tStep, wave.GetData(i));
    }

    // delete f;
  } // end for

  if (debugMode)
  {
    // ----------------------------------------------------------
    // 6. Draw
    // ----------------------------------------------------------
    c->Clear();
    c->Divide(1, 2);
    HiMusicWaveForm waveSum;
    int wavesSize = waves.size();
    for (int i = 0; i < wavesSize; i++)
    {
      HiMusicWaveForm w = waves[i];
      waveSum.Add(w);
    }
    // waveSum.Plot(2,"same");
    for (int i = 0; i < nPoints; i++)
    {
      gSum->SetPoint(i, i * tStep, waveSum.GetData(i));
    }
    c->cd(1);
    gOrig->Draw("AL");
    gSum->Draw("L SAME");
    TLegend *leg = new TLegend(0.6, 0.8, 0.9, 0.9);
    leg->AddEntry(gOrig, "Original waveform", "l");
    leg->AddEntry(gSum, "Sum waveform", "l");
    leg->Draw();
    c->cd(2);
    gResidul->Draw("AL");

    // ----------------------------------------------------------
    // 7. Print results
    // ----------------------------------------------------------
    printf("\n----- Found %zu peaks -----\n", amps.size());
    for (size_t i = 0; i < amps.size(); ++i)
      printf("Peak %2zu : amplitude = %8.1f   start-time = %8.1f ns, peak-time:%8.1f ns\n",
             i + 1, amps[i], starts[i], pt[i]);
  }
  return fitwaves;
}

void HiMusicFittingTask::FindPeaks(HiMusicWaveForm &wave, std::vector<int> &peakPositions, int &nPeaks, float &baseline, float &noise_level)
{
  const int nPoints = wave.GetNumOfPoints();
  const float tStep = wave.GetTStep();
  HiMusicWaveForm wave2(nPoints, tStep);
  if (wave.GetType() == 1)
  { /// fType = 1: Cathode signal, need to invert the waveform
    wave.Inverse();
  }
  const float *data = wave.GetData();

  // 1. 鲁棒基线估计（使用中位数和迭代排除信号区域）
  std::vector<float> dataVec(data, data + nPoints);
  sort(dataVec.begin(), dataVec.end());
  baseline = dataVec[nPoints / 2];

  // 迭代优化基线估计
  float prevBaseline = 0;
  int iteration = 0;
  while (fabs(baseline - prevBaseline) > 1.0 && iteration < 10)
  {
    prevBaseline = baseline;
    std::vector<float> baselinePoints;
    for (int i = 0; i < nPoints; i++)
    {
      // 排除信号区域（超过基线+3倍噪声）
      if (data[i] < baseline + 10)
      {
        baselinePoints.push_back(data[i]);
      }
    }
    if (!baselinePoints.empty())
    {
      sort(baselinePoints.begin(), baselinePoints.end());
      baseline = baselinePoints[baselinePoints.size() / 2];
    }
    iteration++;
  }

  // 2. 噪声估计（使用MAD）
  std::vector<float> deviations;
  for (int i = 0; i < nPoints; i++)
  {
    if (data[i] < baseline + 10)
    { // 只考虑非信号区域
      deviations.push_back(fabs(data[i] - baseline));
    }
  }
  sort(deviations.begin(), deviations.end());
  float MAD = deviations.size() > 0 ? deviations[deviations.size() / 2] : 1.0;
  noise_level = MAD / 0.6745;

  // 3. 计算一阶差分（检测上升沿）
  std::vector<float> diff;
  for (int i = 0; i < nPoints - 1; i++)
  {
    diff.push_back(data[i + 1] - data[i]);
  }

  // 4. 平滑差分信号（5点移动平均）
  std::vector<float> smooth_diff;
  if (diff.size() > 4)
  {
    // 前两点
    smooth_diff.push_back(diff[0]);
    smooth_diff.push_back((diff[0] + diff[1] + diff[2]) / 3.0);

    // 中间点
    int diffSize = diff.size();
    for (int i = 2; i < diffSize - 2; i++)
    {
      smooth_diff.push_back((diff[i - 2] + diff[i - 1] + diff[i] + diff[i + 1] + diff[i + 2]) / 5.0);
    }

    // 最后两点
    smooth_diff.push_back((diff[diff.size() - 3] + diff[diff.size() - 2] + diff[diff.size() - 1]) / 3.0);
    smooth_diff.push_back(diff[diff.size() - 1]);
  }
  else if (!diff.empty())
  {
    smooth_diff = diff; // 数据点太少，不做平滑
  }

  // 4.5 估计平滑差分信号的噪声水平（新增）
  float diffNoiseLevel = 1.0; // 默认值
  if (!smooth_diff.empty())
  {
    std::vector<float> diffNoisePoints;
    float signalThreshold = 3 * noise_level; // 使用原始噪声定义非信号区域

    int smoothDiffSize = smooth_diff.size();
    for (int i = 0; i < smoothDiffSize; i++)
    {
      // 只使用原始波形中非信号区域的点（基线附近）
      if (fabs(data[i] - baseline) < signalThreshold &&
          fabs(data[i + 1] - baseline) < signalThreshold)
      {
        diffNoisePoints.push_back(smooth_diff[i]);
      }
    }

    // 如果非信号区域点太少，使用所有点
    if (diffNoisePoints.size() < 100)
    {
      diffNoisePoints.assign(smooth_diff.begin(), smooth_diff.end());
    }

    // 计算MAD估计噪声
    std::vector<float> diffDeviations;
    for (float val : diffNoisePoints)
    {
      diffDeviations.push_back(fabs(val));
    }
    sort(diffDeviations.begin(), diffDeviations.end());
    float diffMAD = diffDeviations[diffDeviations.size() / 2];
    diffNoiseLevel = diffMAD / 0.6745;
  }

  // 5. 设置检测阈值
  const float thr_rise = 3.0 * noise_level; // 上升沿阈值
                                            // const float thr_rise = 10.0 * diffNoiseLevel;  // 使用差分噪声水平
                                            // const float thr_rise = 3.0 * diffNoiseLevel;  // 使用差分噪声水平
  const int min_rise_points = 4;            // 最小连续上升点数 (100ns)
  const int min_peak_distance = 10;         // 最小峰间距(250ns)
  const int search_window = 20;             // 峰值搜索窗口(250ns)

  // 6. 检测上升沿
  std::vector<int> riseEdges;
  int rise_count = 0;
  int ssmoothDiffSize = smooth_diff.size();
  for (int i = 0; i < ssmoothDiffSize; i++)
  {
    if (smooth_diff[i] > thr_rise)
    {
      rise_count++;
      if (rise_count >= min_rise_points)
      {
        // 记录上升沿起始点
        if (rise_count == min_rise_points)
        {
          riseEdges.push_back(i - min_rise_points + 1);
        }
      }
    }
    else
    {
      rise_count = 0;
    }
  }

  // 7. 在上升沿后寻找峰值
  std::vector<int> candidatePeaks;
  for (int edge : riseEdges)
  {
    int start = std::max(0, edge);
    int end = std::min(nPoints - 1, edge + search_window);

    float maxVal = data[start];
    int maxPos = start;
    for (int j = start; j <= end; j++)
    {
      if (data[j] > maxVal)
      {
        maxVal = data[j];
        maxPos = j;
      }
    }
    candidatePeaks.push_back(maxPos);
  }
  if (debugMode)
  {
    LOG(debug) << "candidatePeaks.size(): " << candidatePeaks.size();
  }
  // 8. 去除重复峰
  sort(candidatePeaks.begin(), candidatePeaks.end());
  std::vector<int> uniquePeaks;
  if (!candidatePeaks.empty())
  {
    uniquePeaks.push_back(candidatePeaks[0]);
    int candidatePeaksSize = candidatePeaks.size();
    for (int i = 1; i < candidatePeaksSize; i++)
    {
      if (candidatePeaks[i] - uniquePeaks.back() > min_peak_distance)
      {
        uniquePeaks.push_back(candidatePeaks[i]);
      }
    }
  }
  if (debugMode)
  {
    LOG(debug) << "uniquePeaks.size(): " << uniquePeaks.size();
  }
  // 9. 基于形状分析区分相邻峰
  std::vector<int> finalPeaks;
  int uniquePeaksSize = uniquePeaks.size();
  for (int i = 0; i < uniquePeaksSize; i++)
  {
    int peakPos = uniquePeaks[i];

    // 检查是否为孤立峰
    if (finalPeaks.empty())
    {
      finalPeaks.push_back(peakPos);
      continue;
    }

    // 检查与前一个峰的距离
    int prevPeak = finalPeaks.back();
    int distance = peakPos - prevPeak;

    // 如果距离足够大，直接添加
    if (distance > 40)
    { // 1000ns
      finalPeaks.push_back(peakPos);
      continue;
    }

    // 分析两个峰之间的谷底
    float valleyValue = data[prevPeak];
    for (int j = prevPeak + 1; j < peakPos; j++)
    {
      if (data[j] < valleyValue)
      {
        valleyValue = data[j];
      }
    }

    // 计算峰谷深度比
    float peak1Height = data[prevPeak] - valleyValue;
    float peak2Height = data[peakPos] - valleyValue;
    float minPeakHeight = std::min(peak1Height, peak2Height);

    // 计算谷底深度（相对于基线）
    // float valleyDepth = baseline - valleyValue;

    // 区分相邻峰的条件
    if (minPeakHeight > 10 * noise_level && // 两个峰都有足够高度
                                            //            valleyDepth > 5 * noise_level &&     // 谷底足够深
                                            // minPeakHeight > 0.3 * max(peak1Height, peak2Height) && // 小峰高度足够
        distance > 4)
    { // 最小距离100ns
      // 确认是两个独立的峰
      finalPeaks.push_back(peakPos);
    }
    else
    {
      // 可能是一个峰或主峰，保留较高的峰
      if (data[peakPos] > data[prevPeak])
      {
        finalPeaks.pop_back();
        finalPeaks.push_back(peakPos);
      }
    }
  } // end for
  if (debugMode)
  {
    LOG(debug) << "finalPeaks.size(): " << finalPeaks.size();
  }

  // 10. 过滤低幅度假峰
  std::vector<int> filteredPeaks;
  for (int pos : finalPeaks)
  {
    if (data[pos] - baseline > 5 * noise_level)
    {
      filteredPeaks.push_back(pos);
    }
  }

  // 11. 设置输出
  nPeaks = filteredPeaks.size();
  peakPositions = filteredPeaks;

  // 调试信息
  if (debugMode)
  {
    LOG(debug) << "===== Advanced Peak Finding Results =====";
    LOG(debug) << "Baseline: " << baseline << " (estimated with " << iteration << " iterations)";
    LOG(debug) << "Noise level: " << noise_level << ", Difference noise level: " << diffNoiseLevel;
    LOG(debug) << "Rise threshold: " << thr_rise;
    LOG(debug) << "Initial candidates: " << candidatePeaks.size();
    LOG(debug) << "Final peaks: " << nPeaks;
    for (int i = 0; i < nPeaks; i++)
    {
      float peakTime = peakPositions[i] * tStep;
      LOG(debug) << "Peak " << i + 1 << ": Position = " << peakPositions[i]
                 << ", Time = " << peakTime << " ns, Amplitude = "
                 << data[peakPositions[i]] - baseline;
    }

    // 绘图
    TCanvas *c1 = new TCanvas("c1", "Advanced Peak Detection", 1200, 800);
    c1->Divide(1, 2);

    // 原始波形
    c1->cd(1);
    TGraph *grWave = new TGraph(nPoints);
    for (int i = 0; i < nPoints; i++)
    {
      grWave->SetPoint(i, i * tStep, data[i]);
    }
    grWave->SetTitle("Waveform with Detected Peaks");
    grWave->GetXaxis()->SetTitle("Time (ns)");
    grWave->GetYaxis()->SetTitle("Amplitude");
    grWave->Draw("AL");

    // 标记峰位
    TGraph *peaks = new TGraph();
    for (int i = 0; i < nPeaks; i++)
    {
      peaks->SetPoint(i, peakPositions[i] * tStep, data[peakPositions[i]]);
    }
    peaks->SetMarkerStyle(23);
    peaks->SetMarkerColor(2);
    peaks->SetMarkerSize(1.);
    peaks->Draw("P SAME");

    // 标记基线
    TLine *baselineLine = new TLine(0, baseline, nPoints * tStep, baseline);
    baselineLine->SetLineColor(8);
    baselineLine->SetLineStyle(2);
    baselineLine->Draw();

    // 差分信号
    c1->cd(2);
    TGraph *grDiff = new TGraph(smooth_diff.size());
    int sssmoothDiffSize = smooth_diff.size();
    for (int i = 0; i < sssmoothDiffSize; i++)
    {
      grDiff->SetPoint(i, i * tStep, smooth_diff[i]);
    }
    grDiff->SetTitle("Smoothed Difference");
    grDiff->GetXaxis()->SetTitle("Time (ns)");
    grDiff->GetYaxis()->SetTitle("Difference");
    grDiff->Draw("AL");

    // 阈值线
    TLine *thresholdLine = new TLine(0, thr_rise, nPoints * tStep, thr_rise);
    thresholdLine->SetLineColor(2);
    thresholdLine->SetLineStyle(2);
    thresholdLine->Draw();

    // 标记上升沿
    TGraph *riseMarkers = new TGraph();
    int riseEdgesSize = riseEdges.size();
    for (int i = 0; i < riseEdgesSize; i++)
    {
      riseMarkers->SetPoint(i, riseEdges[i] * tStep, thr_rise * 1.1);
    }
    riseMarkers->SetMarkerStyle(22);
    riseMarkers->SetMarkerColor(4);
    riseMarkers->SetMarkerSize(1.);
    riseMarkers->Draw("P SAME");
  }
}

HiMusicWaveForm HiMusicFittingTask::CreateWave(TF1 *f, int np, int dt)
{
  HiMusicWaveForm awave;
  awave.SetNPoints(np);
  awave.SetTStep(dt);

  // int Npar = f->GetNpar();
  double *par = f->GetParameters();
  TString fname = f->GetName();
  if (fname != "FitSingleExp" && fname != "FitBiExp" && fname != "FitLineExp")
  {
    LOG(error) << "fit type error !!! Only support singleExp, biExp and lineExp.";
    return awave;
  }
  double data;
  double x[1];
  for (int i = 0; i < np; i++)
  {
    x[0] = i * dt;
    if (fname == "FitSingleExp") /// single Exp
    {
      data = HiMusicFittingTask::singleExp(x, par);
    }
    else if (fname == "FitBiExp") /// bi Exp
    {
      if (x[0] < f->GetX(0))
        data = 0.0;
      else
        data = HiMusicFittingTask::biExp(x, par);
    }
    else if (fname == "FitLineExp")
    {
      data = HiMusicFittingTask::lineExp(x, par);
    }
    awave.SetData(i, data);
  }
  return awave;
}

double HiMusicFittingTask::singleExp(double *x, double *par) ////single Exponential with pedestal substracted
{
  /// par[0]: peak height
  /// par[1]: peak position
  /// par[2]: tao
  return par[0] * exp(-(x[0] - par[1]) / par[2]);
}

double HiMusicFittingTask::biExp(double *x, double *par) ////bi-Exponential with pedestal substracted
{
  // par[0] = A  (amplitude)
  // par[1] = t0 (start time)
  // par[2] = τ1 (rise)
  // par[3] = τ2 (decay)
  double t = x[0];
  double t0 = par[1];
  if (t < t0)
    return 0.0;
  double tau1 = par[2];
  double tau2 = par[3];
  double val = par[0] * (exp(-(t - t0) / tau2) - exp(-(t - t0) / tau1));
  return val;
}

double HiMusicFittingTask::lineExp(double *x, double *par) ////Linear rise + Exponential decay with pedestal substracted
{
  // par[0] = A  (amplitude)
  // par[1] = t0 (start time)
  // par[2] = tp (peak time)
  // par[3] = τ (decay time)
  Double_t retval = 0;
  if (x[0] > par[1] && x[0] < par[2])
    retval = par[0] * (x[0] - par[1]) / (par[2] - par[1]);
  else if (x[0] >= par[2])
    retval = par[0] * exp(-(x[0] - par[2]) / par[3]);
  return retval;
}